#####################################################################
# Figure 4, Tables 4-5
#####################################################################

#####################################################################
# Preliminaries
#####################################################################

rm(list = ls())

library(xtable)
library(limSolve)
library(synthdid)
library(foreign)
library(scinference)
library(writexl)

set.seed(12345)

source("code/auxiliary scripts/common_functions.R")

###################################################################
# Data (data source: replication package of Andersson (2019))
###################################################################

# Data prep

carbontax_data <- read.dta("data/carbontax_data.dta")

keeps <- c("Countryno", "year","CO2_transport_capita")
carbontax_data <- carbontax_data[keeps]

carbontax_wide <- reshape(carbontax_data, timevar = "year",idvar = c("Countryno"),direction = "wide")

Y1 <- as.matrix(cbind(carbontax_wide[13,-1]))
Y0 <- t(as.matrix(carbontax_wide[-13,-1]))

T0  <- 30 # 1960-1989
T1  <- length(Y1)-T0
T01 <- T0 + T1

J <- dim(Y0)[2]

# Figure 4

graphics.off()
time <- c(seq(1960,2005,1))
pdf("graphics/sweden_data_raw.pdf",pointsize=14,width=9.0,height=6.0)
plot(NULL, ylab="Per capita CO2 emissions", xlab="Year", main="", col="blue", pch=".", xlim = range(time), ylim=c(0,7))
for (j in 1:dim(Y0)[2]) lines(time,Y0[,j],col="darkgray",lwd=0.5,lty=1)
lines(time,Y1,col="black",lwd=5,lty=1)
legend("topleft",legend=c("Sweden","Controls"), seg.len=2, col=c("black","darkgray"),fill=NA,border=NA, lty=c(1,1),lwd=c(5,0.5), merge=T,bty="n")
dev.off()

###################################################################
# Placebo test
###################################################################

alpha.sig <- 0.1

Y1.pre <- Y1[1:T0]
Y0.pre <- Y0[(1:T0),]

# Placebo test with K=3

r.placebo.1812  <- scinference(Y1.pre,Y0.pre,T1=12,T0=18,inference_method="ttest",estimation_method="sc",K=3,alpha=alpha.sig,lsei_type = 1)
r.placebo.219   <- scinference(Y1.pre,Y0.pre,T1=9,T0=21,inference_method="ttest",estimation_method="sc",K=3,alpha=alpha.sig,lsei_type = 1)

# Table 4

tau.placebo <- c(r.placebo.1812$att,r.placebo.219$att)
lb.placebo  <- c(r.placebo.1812$lb,r.placebo.219$lb)
ub.placebo  <- c(r.placebo.1812$ub,r.placebo.219$ub)

print("LaTeX output corresponding to Table 4")
results_table4 <- cbind(tau.placebo,lb.placebo,ub.placebo)
print(xtable(cbind(c(18,21),results_table4),digits=c(0,0,2,2,2)),include.rownames = F)

data_results_table4 <- as.data.frame(matrix(sprintf("%.2f",results_table4),2,3))
colnames(data_results_table4) <- NULL
write_xlsx(data_results_table4,"tables/table4.xlsx")

###################################################################
# Main results
###################################################################

alpha.sig <- 0.1

# t-Test with synthetic control, K=3
r.sc.3    <- scinference(Y1,Y0,T1=T1,T0=T0,inference_method="ttest",estimation_method="sc",K=3,alpha=alpha.sig,lsei_type = 1) 
tau.sc.3  <- r.sc.3$att
se.sc.3   <- r.sc.3$se
lb.sc.3   <- r.sc.3$lb  
ub.sc.3   <- r.sc.3$ub  

# t-Test with DiD, K=3
r.did   <- scinference(Y1,Y0,T1=T1,T0=T0,inference_method="ttest",estimation_method="did",K=3,alpha=alpha.sig,lsei_type = 1) 
tau.did <- r.did$att
se.did  <- r.did$se
lb.did  <- r.did$lb  
ub.did  <- r.did$ub  

# Synthetic DiD
tau.sdid  <- synthdid_estimate(t(cbind(Y0,c(Y1))),J,T0)
se.sdid   <- sqrt(vcov(tau.sdid, method='placebo'))
lb.sdid   <- tau.sdid - qnorm(1-alpha.sig/2)*se.sdid 
ub.sdid   <- tau.sdid + qnorm(1-alpha.sig/2)*se.sdid 

# Li (2020)'s subsampling approach
r.li    <- li2020(Y1,Y0,T0=T0,T1=T1,20,1000,alpha.sig)
tau.li  <- r.li$ATT_SC
lb.li   <- r.li$ci_low
ub.li   <- r.li$ci_high

# Table 5

print("LaTeX output corresponding to Table 5")

results_table5 <- c(tau.sc.3,lb.sc.3,ub.sc.3,tau.did,lb.did,ub.did,tau.li,lb.li,ub.li,tau.sdid,lb.sdid,ub.sdid)
print(xtable(rbind(results_table5),digits=c(0,rep(2,12))),include.rownames = F)

data_results_table5 <- as.data.frame(sprintf("%.2f",results_table5))
colnames(data_results_table5) <- NULL
write_xlsx(data_results_table5,"tables/table5.xlsx")

